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1.  Introduction 


Nitromethane  (NM)  is  one  of  the  most  widely  studied  insensitive  energetic 
molecules.1-8  The  main  reason  for  its  popularity  is  that  it  is  one  of  the  simplest 
nitro-organic  molecules  containing  only  one  nitro  group  (NO2)  bonded  to  one 
methyl  group  (CH3)  through  a  C-N  bond.  This  makes  it  an  ideal  candidate  to  study 
chemical  reaction  mechanisms  associated  with  conventional  explosive  initiation 
and  subsequent  detonation.3,6,7  However,  understanding  the  ultra-fast  reaction 
mechanisms  for  initiation  is  challenging. 

The  study  of  initiation  and  chemical  reaction  mechanisms  of  mechanically  stressed 
materials  such  as  shock  compression,  hydrostatic  compression,  uniaxial 
compression,  or  shear  is  called  mechanochemistry.  The  broad  field  of 
mechanochemistry  applied  to  energetics  covers  a  large  range  of  conditions  and 
interactions  that  can  be  potentially  exploited  for  advanced  energy  release 
applications.  In  addition,  mechanochemistry  is  also  relevant  for  nonenergetics  for 
phase  transitions  and  chemical  reactions  that  occur  when  mechanically  stressed. 

NM  is  a  liquid  at  ambient  pressure,  but  when  compressed  to  about  0.3  GPa,  it 
crystallizes  in  the  P2i2i2i  space  group.1  When  hydrostatically  compressed  to 
28  GPa,  NM  slowly  decomposes  to  what  is  reported  to  be  a  3-D  disordered 
polymer.2  In  contrast,  when  NM  is  compressed  to  28  GPa  and  then  sheared,  Raman 
studies  show  that  the  decomposition  is  sudden  and  explosive.5  In  addition, 
structural  modifications  are  found  to  occur  at  a  lower  hydrostatic  pressure  when 
sheared  in  comparison  to  hydrostatic  pressure  alone.5  One  model  that  has  been 
proposed  to  explain  the  explosive  decomposition  of  NM  under  shear  stress 
examines  the  orientation-dependent  steric  hindrance  to  shear  flow  at  the  molecular 
level.4  It  was  found  that  certain  orientations  of  NM  are  highly  sterically  hindered 
and  are  therefore  more  sensitive  to  chemical  reactions  and  initiation  than  others. 
However,  this  is  a  static  model  that  does  not  capture  the  dynamics  that  occur  in  the 
crystal  during  shear.4 

In  this  work  we  examine  the  dynamics  that  occur  during  shear  of  NM  by  performing 
classical  molecular  dynamics  simulations.  The  crystallographic  orientation 
dependence  on  the  shear  response  of  NM  is  calculated.  In  addition,  since  there  is  a 
strong  correlation  between  grain  boundaries  (GBs)  and  chemical  initiation  in 
energetic  materials,  bicrystalline  NM  with  a  GB  is  also  simulated. 
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2.  Computational  Details 


Molecular  dynamics  simulations  are  perfonned  using  the  LAMMPS  (Large-scale 
Atomic/Molecular  Massively  Parallel  Simulator)  simulation  package.9  The 
interatomic  interactions  are  described  by  the  ReaxFF-lg  potential  recently  created 
using  the  MOES  (Multiple  Objective  Evolutionary  Strategies)  algorithm1011  and 
refitted  to  NM.10  13  This  potential  shows  good  agreement  to  the  experimental 
equation  of  state  and  melting  lines  for  NM.  Fully  3-D  periodic  crystals  of  NM  in 
various  orientations  are  generated  using  the  GB  builder  code  developed  by  Shawn 
Coleman  (US  Army  Research  Laboratory,  Weapons  and  Materials  Research 
Directorate)  and  initially  relaxed  in  the  NPT  (constant  number  of  atoms,  pressure, 
and  temperature)  ensemble  at  20  GPa  and  298  K.  Special  care  is  taken  to  ensure  the 
crystal  structure  generated  by  the  code  is  equivalent  to  the  original  crystal  structure. 
The  crystals  are  then  sheared  in  the  microcanonical  ensemble  with  a  constant  strain 
rate  equal  to  0.0008  ps"1  over  a  100-ps  time  period.  The  strain  rate  is  defined  as 
offset  divided  by  the  length,  where  the  length  is  the  initial  box  length  in  the  direction 
perpendicular  to  shear  (e.g.,  the  y  box  length  for  xy  deformation),  and  the  offset  is 
the  change  in  the  direction  in  the  shear  direction  (e.g.,  the  change  in  the  x  box  length 
for  xy  deformation).  Each  crystal  is  sheared  along  all  3  angles  in  the  simulation  box 
(xy,  xz,  and  yz).  The  effect  the  strain  rate  has  on  the  results  is  analyzed  by 
performing  simulations  with  a  strain  rate  an  order  of  magnitude  higher  and  lower 
than  0.0008,  and  no  significant  deviations  are  found.  Each  orientation  has  a 
different  number  of  atoms  in  the  unit  cell  to  ensure  3-D  periodicity;  however,  all 
simulations  use  at  least  20,000  atoms  in  a  simulation  box  bigger  than  35  x  35  x 
35  A3.  The  number  of  atoms  in  the  cell  is  also  varied  to  ensure  no  significant  size 
dependence  on  the  results.  Several  hundred  bicrystalline  structures  of  NM  are 
constructed  using  similar  orientations  used  in  the  single  crystal  simulations  and  are 
initially  relaxed  at  20  GPa  and  298  K.  The  bicrystalline  structures  are  2-D  periodic 
in  the  xz  plane  while  the  size  of  the  grain  along  the  y  axis  is  100  A  (both  grains  are 
100  A).  Simulations  with  50  A  and  300  A  are  found  to  give  similar  results. 

3.  Results 


A  schematic  for  the  shear  simulations  is  shown  in  Fig.  la  and  lb  for  orientation 
(Oil)  deformed  along  the  xy  angle.  The  (Oil)  orientation  is  directed  along  the  z 
axis.  Since  NM  has  an  orthorhombic  unit  cell,  each  cell  angle  is  initially  90°.  As 
the  crystal  is  sheared,  the  xy  angle  decreases  while  the  atoms  move  simultaneously 
according  to  Newton’s  laws  to  minimize  the  energy  and  forces. 
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Fig.  1  Schematic  of  the  shear  simulations  on  NM  using  the  single  crystal  (Oil)  orientation 
sheared  along  the  xy  angle  as  an  example.  Initially  a)  each  angle  of  the  simulation  box  is  90° 
and  then  b)  the  xy  angle  is  reduced  to  simulate  shearing  the  material.  In  b)  NM  has 
transformed  into  a  disordered  state. 

The  response  of  NM  under  shear  is  highly  dependent  on  the  crystallographic 
orientation.  In  Fig.  2a  the  shear  stress  is  given  as  a  function  of  shear  angle  for  the 
(Oil)  orientation  for  the  3  shear  directions.  This  orientation  is  representative  of  the 
results  found  for  many  other  orientations  that  are  simulated.  The  shear  stress  profile 
is  different  for  all  3  shear  directions,  but  in  all  cases  the  shear  stress  builds  up  until 
a  maximum  shear  stress  is  reached.  After  the  maximum  shear  stress  is  reached,  the 
material  suddenly  transforms  into  a  disordered  state  as  shown  in  Fig.  lb  and  the 
shear  stress  suddenly  drops  (see  Fig.  2a  for  shear  along  the  xy  angle).  In  addition 
to  the  drop  in  the  shear  stress,  the  transformation  to  the  disordered  state  is 
accompanied  with  a  drop  in  potential  energy  as  shown  in  Fig.  2b,  which  is  most 
easily  seen  for  the  shear  of  the  xy  angle  at  about  65°.  The  drop  in  potential  energy 
corresponds  to  a  sudden  increase  in  kinetic  energy.  Nevertheless,  the  rapid  increase 
in  kinetic  energy  does  not  cause  reactions  to  occur.  The  magnitude  of  the  maximum 
shear  stress,  the  shear  angle,  and  the  magnitude  of  the  potential  energy  drop  are  all 
highly  dependent  on  the  crystallographic  orientation  and  shear  direction. 
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Fig.  2  The  a)  shear  stress  Si  (xy,  xz,  or  yz)  and  b)  potential  energy  as  a  function  of  shear  angle 
(xy,  xz,  or  yz)  during  the  shear  simulation  for  the  (Oil)  orientation  of  NM.  This  orientation 
is  chosen  to  be  displayed  as  a  representative  example  of  a  large  shear  stress  and  large  shear 
angle  as  well  as  a  low  shear  stress  and  low  shear  angle  that  occurs  in  NM  during  shear. 

The  orientation  dependence  can  be  attributed  to  the  atomistic  dynamics  that  occur 
in  the  crystal  during  shear.  For  shear  along  the  xy  angle,  the  NM  molecules  are 
initially  oriented  so  that  they  are  staggered  with  respect  to  one  another  along  the 
y-axis  as  shown  in  Fig.  3a.  As  the  crystal  is  sheared,  the  NM  molecules  rotate  so 
that  they  are  roughly  oriented  in  the  same  direction,  as  can  be  seen  in  Fig.  3b  at 
65.77°,  which  is  just  before  the  transformation  to  the  disordered  state.  The  rotation 
of  the  molecules  is  what  allows  the  crystal  to  be  sheared  to  such  a  large  angle,  as 
shown  in  Fig.  2a.  Once  the  molecules  can  no  longer  rotate,  a  large  shear  stress  of 
about  4  GPa  builds  up,  and  then  the  transformation  to  the  disordered  state  occurs. 

(a)  xy  z.  =90°  (b)  xy  z  =65.77° 


Fig.  3  Zoomed-in  snapshots  of  the  NM  molecules  for  the  (Oil)  orientation  sheared  along  the 
xy  angle.  The  NM  molecules  are  a)  initially  staggered  along  the  y  axis,  and  then  after  the 
crystal  is  sheared  b),  the  NM  molecules  are  rotated  to  orient  themselves  in  the  same  direction. 

In  contrast,  when  the  (Oil)  orientation  is  sheared  along  the  yz  angle,  only  a  small 
shear  stress  is  built  up  (about  2  GPa),  as  shown  in  Fig.  2a.  This  is  because  in 
addition  to  the  rotation  of  the  NM  molecules,  the  planes  of  the  NM  molecules  also 
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shift  in  the  shear  direction,  reducing  the  shear  stress.  The  initial  configuration  of 
the  NM  molecules  is  displayed  in  Fig.  4a,  where  the  NM  molecules  are  once  again 
staggered  along  the  y  axis.  Notice  that  there  is  some  empty  space  between  the 
“planes”  of  NM  molecules  in  the  xy  plane  in  Fig.  4a.  This  allows  the  top  plane  of 
NM  molecules  to  shift  with  respect  to  the  bottom  plane  during  shear.  This  can  be 
seen  by  comparing  the  relative  positions  of  the  NM  molecules  in  Fig  4a  with  those 
in  Fig.  4b,  where  the  y-coordinate  of  adjacent  NM  molecules  along  the  z  axis  are 
no  longer  equivalent. 


(a)  yz  L  =90°  (b)  yz  z  =72.77° 


Fig.  4  Zoomed-in  snapshots  of  the  NM  molecules  for  the  (Oil)  orientation  sheared  along  the 
yz  angle.  The  NM  molecules  are  initially  a)  staggered  along  the  y-axis,  and  the  adjacent  planes 
of  NM  molecules  have  approximately  the  same  y  coordinate.  After  shear  b),  the  NM  molecules 
are  oriented  along  the  same  direction,  and  the  adjacent  plane  of  NM  has  shifted  along  the  y 
axis. 

For  all  the  orientations  simulated,  the  maximum  shear  stress  ranges  between  1.68 
and  5.24  GPa,  and  the  shear  angle  ranges  from  64.3°  to  86.7°  as  shown  in  Table  1. 
The  (Oil)  orientation  given  as  an  example  has  a  range  of  2.35-4.48  GPa  and 
65.5°-81.5°  for  the  3  shear  directions,  so  it  is  a  good  representative  example  of  the 
variety  of  shear  dynamics  in  NM.  The  orientations  that  can  be  sheared  to  the  lowest 
angles  likely  have  a  significant  amount  of  rotation  in  the  NM  molecules,  such  as 
what  is  observed  in  the  (Oil)  orientation.  In  contrast,  the  orientations  that  reach  the 
maximum  stress  near  90°  likely  do  not  allow  for  any  rotation  of  the  molecules.  At 
this  time  it  is  not  clear  whether  initiation  is  preferred  for  orientations  with  the 
highest/lowest  shear  stress  or  highest  shear  angle,  or  if  another  parameter  such  as 
maximum  temperature  controls  the  sensitivity  in  each  orientation. 
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Table  1  Summary  of  the  maximum  shear  stress  and  the  shear  angle  at  maximum  shear  stress 
for  each  single  crystal  orientation  of  NM  simulated 


Orientation 

(hkl) 

Max  |S„( 
(GPa) 

Shear  angle 

o 

Max  |Sxz| 
(GPa) 

Shear  angle 

o 

Max  |Syz| 
(GPa) 

Shear 

angle 

(°) 

(Oil) 

4.21 

83.1 

5.24 

80.8 

4.03 

81.5 

(101) 

2.15 

69.1 

2.70 

85.2 

4.67 

82.1 

(111) 

3.73 

83.3 

3.34 

82.9 

4.54 

82.1 

(110) 

3.03 

66.9 

4.53 

80.8 

2.98 

85.2 

(010) 

3.92 

64.7 

3.46 

85.1 

1.95 

72.9 

(111) 

4.06 

82.6 

3.03 

86.7 

2.96 

77.1 

(111) 

4.08 

82.7 

3.27 

76.1 

3.01 

86.7 

(110) 

3.43 

68.4 

3.46 

86.1 

2.85 

85.1 

(Oil) 

4.38 

65.5 

4.48 

81.5 

2.35 

71.6 

(101) 

2.19 

64.3 

2.79 

84.9 

3.51 

78.3 

(Oil) 

3.07 

67.4 

2.94 

86.6 

1.68 

72.9 

(001) 

2.60 

67.5 

2.77 

72.7 

4.27 

64.9 

(100) 

4.26 

64.6 

3.59 

66.2 

2.54 

64.2 

In  order  to  include  GBs  into  the  simulations,  the  structure  at  the  GB  interface  should 
be  similar  to  what  is  present  in  experiment;  however  the  experimental  GB  structure 
of  NM  is  currently  unknown.  To  generate  realistic  GB  structures,  the  following 
method  is  used.  Several  hundred  GB  structures  are  generated  by  creating  3  GBs 
with  different  misorientation  angles  (along  z  axis)  for  grains  with  orientations 
(h1/c1i1)/(/i2/c2^2)  (along  the  y  axis)  where  h,  k  and  /  range  from  -1  to  1.  A  total  of 
273  GBs  are  created;  however,  some  of  them  have  the  same  (hkl)  orientation  and 
therefore  represent  a  bulk  structure  with  a  defect  at  the  interface  rather  than  a  real 
GB.  Nevertheless,  all  the  structures  are  relaxed  at  20  GPa  and  300  K  and  ranked 
according  to  the  GB  energy.  The  following  fonnula  is  used  for  the  GB  energy: 
(UCb  ~  U perf  ect) /2A,  where  Ugh  is  the  GB  energy  and  Uperfect  is  the  energy  of  a 
perfect,  single  crystal,  and  A  is  the  area  of  the  GB.  The  shear  simulations  are 
performed  on  the  1 5  best  structures  with  the  lowest  GB  energy. 

When  a  bicrystalline  GB  is  included  in  the  shear  simulations,  the  shear  dynamics 
are  more  complex;  however,  the  overall  process  is  still  similar  to  the  single  crystal. 
The  shear  stress  still  builds  up  followed  by  a  transfonnation  to  a  disordered  state. 
In  Fig.  5  a  and  5b,  the  shear  stress  and  potential  energy  profile  during  shear  are  given 
as  a  function  of  shear  angle  for  the  (011)/(011)  orientation  along  the  z  axis  and 
(0  27  14)/(0  27  14)  orientation  along  the  y  axis.  This  corresponds  to  a 
misorientation  angle  of  90°  about  the  z  axis.  This  orientation  is  chosen  as  an 
example  because  it  has  similar  properties  as  many  of  the  GBs  that  are  simulated. 
As  shown  in  Table  2,  the  maximum  shear  stresses  and  the  GB  energy  for  this 
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orientation  are  similar  to  many  of  the  orientations  simulated,  aside  from  a  few 
outliers. 


Angle  (degrees)  Angle  (degrees) 


Fig.  5  The  (a)  shear  stress  and  (b)  potential  energy  given  as  a  function  of  shear  angle  for  a 
bicrystalline  simulation  of  NM 


Table  2  Summary  of  the  calculated  grain  boundary  (GB)  orientations,  maximum  shear 
stress,  shear  angle  where  the  shear  stress  occurs,  and  GB  energy 


Z-axis  orientation: 

( hkl)/(hkl ) 

Y-axis  orientation: 
(hkl)/(hkl) 

GB 

energy 

(eV/A2) 

Max. 

|Sxy| 

(GPa) 

Angle 

(°) 

Max. 

|Sxz| 

(GPa) 

Angle 

(°) 

Max. 

|Sxz| 

(GPa) 

Angle 

(°) 

<011>/<011> 

<0  27  14>/<0  27  14) 

0.903 

1.56 

85.0 

1.91 

76.6 

2.10 

70.6 

<011)/<011) 

<0  27  14)/<0  27  14) 

1.025 

1.46 

70.6 

1.76 

85.7 

2.11 

71.6 

<011)/<011) 

<0  27  14)/<0  27  14) 

0.932 

1.33 

79.5 

1.70 

82.5 

2.87 

60.9 

<110)/<110) 

<3  2  0)/(49  33  0) 

0.867 

1.80 

85.3 

3.64 

73.8 

1.64 

59.7 

<111>/<111) 

<8  15  5>/<8  15  5) 

0.819 

2.55 

77.2 

2.48 

82.2 

2.06 

68.7 

<111)/<111) 

<8  15  5)/(8  15  5) 

1.246 

2.09 

78.1 

4.15 

83.1 

1.92 

69.7 

<111>/<111) 

<8  15  5)/<8  15  5) 

3.111 

2.09 

80.9 

2.55 

64.2 

1.84 

82.9 

<111>/<111) 

<8  15  5)/<8  15  5) 

0.921 

2.14 

77.5 

2.21 

76.6 

2.21 

65.0 

<111)/<111> 

<8  15  5)/(8  15  5) 

0.698 

1.98 

82.6 

2.60 

82.0 

1.98 

69.3 

<111)/<111) 

<8  15  5)/<8  15  5) 

0.987 

1.71 

85.7 

3.76 

83.4 

2.17 

67.5 

<111>/<111> 

<8  15  5)/<8  15  5) 

0.773 

2.00 

81.6 

2.34 

74.9 

1.92 

69.5 

<111)/<111) 

<8  15  5)/<8  15  5) 

0.770 

1.99 

81.8 

2.62 

81.8 

1.90 

67.5 

<111>/<111> 

<8  15  5)/<8  15  5) 

0.928 

1.85 

83.9 

2.29 

73.5 

1.90 

71.0 

<111)/<111) 

<8  15  5)/<8  15  5) 

0.809 

2.39 

77.9 

2.42 

70.8 

1.91 

68.1 

<011>/<011) 

<0  27  14)/(0  27  14) 

0.449 

1.49 

85.0 

1.73 

72.6 

2.5 

58.4 

As  shown  in  Fig.  5a,  as  this  orientation  is  sheared  along  the  yz  angle,  the  shear 
stress  builds  up  slowly  and  even  drops  slightly  around  85°,  then  builds  up  again.  In 
this  simulation,  the  grain  is  slowly  rotating  with  the  shear,  as  displayed  in  Fig.  6a 
and  6b.  The  molecules  do  not  rotate  much  in  the  beginning  of  the  simulation,  but 
the  grain  rotates  considerably.  This  can  be  seen  more  clearly  by  examining  the 
C-N  bonds  in  NM  along  the  arrow  shown  in  Fig.  6a  and  6b. 
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(a)  Grain  boundary:  <0 1 1  >/<0 1 1  >  along  z.  <0  27  - 1 4>/ <0  -27  14>  along  y 

Angle  89.70325650000001 


25578  atoms 

(b) 

Angle  74.960496 


U: 


25578  atoms 


Fig.  6  Snapshots  during  the  shear  stress  along  the  yz  angle  with  a  bicrystal  GB.  A  line  is 
shown  in  a)  and  b)  as  a  guide  to  the  eye  to  shown  the  rotation  of  the  grain. 


Also  shown  in  Fig.  6a  is  that  the  disordered  material  originates  at  the  GB  interface 
and  expands  from  there  as  shown  in  Fig.  6b.  However,  in  stark  contrast,  when 
sheared  along  the  xz  angle,  the  entire  crystal  suddenly  transforms  into  the 
disordered  state,  which  causes  the  shear  stress  to  drop  rapidly,  as  shown  in  Fig.  5  a 
at  85°.  Interestingly,  the  shear  stress  then  rises  again.  This  is  because  the 
transformation  to  the  disordered  state  is  incomplete,  and  the  material  still  maintains 
a  small  amount  of  order.  The  transformation  to  the  disordered  state  is  not  complete 
until  about  75°. 

As  shown  in  Table  1,  the  largest  absolute  shear  stress  is  found  to  occur  in  the 
(011)/(111)  orientation  along  z  and  (8  15  5)/(8  15  5)  along  y  sheared  along  the 
xz  axis  of  4.15  GPa  at  an  angle  of  83.1°.  However,  in  spite  of  the  large  build-up  of 
shear  stress,  no  reactions  are  observed  in  this  bicrystalline  simulation,  which  is 
similar  to  what  is  found  for  the  single  crystal  simulations. 


4.  Conclusions 


The  short  time-scale  dynamics  during  shear  of  NM  have  been  revealed  using 
molecular  dynamics  simulations.  Shear  simulations  have  been  perfonned  using 
single  crystals  with  a  range  of  orientations  and  shear  directions  as  well  as  bicrystals, 
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also  with  a  range  of  orientations  and  misorientation  angles.  Although  no  chemical 
reactions  are  observed  during  the  simulations,  the  dynamics  presented  are  possibly 
associated  with  the  short  time-scale  reaction  pathway  that  leads  to  initiation.  In  this 
sense,  these  results  are  critical  for  understanding  the  complete  mechanism  for 
initiation  of  NM  under  shear  stress. 

Since  no  reactions  are  observed  in  the  simulations  but  are  observed  in  experiment, 
the  future  work  is  mostly  geared  toward  resolving  this  intriguing  discrepancy  as 
well  as  generally  increasing  the  length  and  time  scales,  and  the  range  of  orientations 
studied.  The  initiation  may  be  related  to  the  presence  of  multiple  GBs,  so 
simulations  with  polycrystalline  NM  may  be  important.  Investigating  the  size- 
dependence  of  the  results  more  thoroughly  with  GBs  may  also  be  needed  because 
the  GB  simulations  require  that  the  size  of  the  grain  be  much  larger  than  the  size  of 
the  GB.  In  addition,  initiation  may  be  related  to  the  slip  and  subsequent  momentum 
transfer  of  grains  much  larger  than  what  is  simulated,  so  the  size-dependence  needs 
to  be  investigated  thoroughly.  The  time  scale  of  the  simulations  may  also  be  an 
issue  because  it  may  be  that  the  time  scale  for  reactions  to  occur  is  simply  larger 
than  the  time  scales  simulated.  An  analysis  of  activation  energies  and  Arrhenius 
reaction  rates  for  proposed  reactions  may  yield  an  estimate  of  the  time  scale  needed. 
A  greater  variety  of  GB  orientations  can  also  be  simulated,  which  may  yield 
interesting  results.  All  the  results  presented  are  at  20  GPa;  it  would  be  interesting 
to  examine  the  pressure  dependence  of  the  results,  since  there  is  a  clear  pressure 
dependence  in  the  experiment  related  to  the  reaction  threshold  pressure. 

This  work  was  conducted  as  part  of  the  DOD  High  Performance  Computing 
Modernization  Program’s  (HPCMP’s)  Workforce  Development  initiative  and  is 
aimed  at  providing  the  future  science  and  engineering  workforce  the  computational 
skills  and  experience  necessary  to  close  the  gap  between  the  technological 
capability  and  computational  skills  necessary  to  support  the  DOD’s  future 
Warfighter  needs.  The  HPCMP’s  2016  Internship  Program  was  a  valuable 
experience  for  Brad  Steele,  and  the  impact  of  the  summer  research  experience  on 
his  professional  career  is  discussed  in  the  Appendix. 
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Appendix.  Impact  of  Summer  Research  Experience 
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The  following  remarks  were  made  by  Brad  Steele  in  reference  to  his  experience  as 
a  student  intern  at  ARL  under  the  2016  HPCMP  Internship  Program. 

Background 

I  am  a  PhD  candidate  at  the  University  of  South  Florida,  Physics  Department,  where 
I  earned  my  MS  in  Applied  Physics.  I  have  a  significant  amount  of  experience 
perfonning  density  functional  theory  calculations  and  bash  scripting  and 
programming  for  running  large  amounts  of  simulations  and  data  analysis.  My  past 
research  experience  dealt  with  predicting  and  characterizing  new  kinds  of  materials 
at  high  pressures.  This  involves  zero-temperature  energy-minimization  calculations 
to  detennine  thermodynamic  and  dynamical  stability.  I  have  much  less  experience 
perfonning  simulations  using  LAMMPS  and  analyzing  time-dependent  dynamical 
results.  My  goals  for  participating  in  the  internship  were  to  gain  more  experience 
using  LAMMPS  and  analyzing  dynamical  results.  Furthennore,  I  wanted  to  gain 
experience  conducting  research  at  a  government  research  lab  to  understand  the 
techniques  and  how  research  challenges  are  tackled  on  a  government  level. 

Impact 

The  HPCMP  Internship  Program  has  been  incredibly  useful  for  me.  I  have  gained 
a  tremendous  amount  of  experience  running  molecular  dynamics  simulations  and 
analyzing  the  results.  I  have  also  been  able  to  further  develop  my  programming  and 
bash  scripting  skills.  More  significantly,  I  have  a  much  better  appreciation  for  what 
molecular  dynamics  is  capable  of  and  what  its  limitations  are.  Before  this 
internship,  I  had  never  used  ReaxFF;  I  now  recognize  its  potential.  Furthermore,  I 
have  a  much  better  understanding  of  the  technical  tenns  used  by  the  energetic 
materials  community  and  what  types  of  properties  they  are  interested  in.  One  thing 
that  surprised  me  about  my  project  was  how  difficult  it  is  to  get  the  metastable 
molecules  to  react.  I  had  generally  assumed  that  chemical  reactions  can  be  observed 
easily  using  molecular  dynamics,  but  that  is  not  the  case.  I  have  also  gained  an 
appreciation  of  what  researchers  mean  when  they  say  a  material  reacts  or 
decomposes  under  some  sort  of  stress,  why  the  response  of  the  material  depends  on 
the  type  of  stress  applied,  and  why  the  rate  of  the  reaction  is  very  important  factor 
for  energetic  materials.  On  a  personal  level,  this  experience  is  of  tremendous  value 
to  me  because  I  want  to  work  at  a  government  lab  when  I  graduate.  It  has  been  very 
good  to  meet  the  other  researchers  and  learn  about  what  they  are  working  on  and 
how  they  are  trying  to  solve  their  research  problems.  This  program  also  fulfills  a 
necessary  requirement  for  my  PhD. 
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List  of  Symbols,  Abbreviations,  and  Acronyms 


GB 

grain  boundary 

HPCMP 

High  Performance  Computing  Modernization  Program 

LAMMPS 

Large-scale  Atomic/Molecular  Massively  Parallel  Simulator 

MOES 

Multiple  Objective  Evolutionary  Strategies 

NM 

nitromethane  (CH2NO3) 
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